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A comprehensive review of all relevant experimental data was completed, including recent data for the drag 
coefficient for a sphere in supersonic and hypersonic flows. The primary characterization parameter included the 
relative Mach, Knudsen, and Reynolds numbers based on the relative velocity, the sphere diameter, and other 
parameters. This review of data showed that the previously proposed nexus at a Reynolds number below 45 was not 
strictly met, and it instead included a weak transonic bump, which was identified numerically for the first time with 
the present simulations. New continuum-gas and rarefied-gas simulations were conducted and were combined with 
the expanded experimental dataset to improve the quantitative description of the drag coefficient in this region. The 
results indicated that a quasi nexus bridges the rarefaction regime and the compressible flow regimes. The 
comprehensive dataset was then used to develop new empirical models for the drag coefficient that showed 
improved robustness and accuracy as compared to previous models. These models are limited by the critical 
Reynolds number associated with boundary-layer transition on the sphere, which was found to increase 


substantially with the sphere Mach number. 


Nomenclature 


= speed of sound in the gas 
= force coefficient 
= velocity coefficient 
= particle diameter 
= force 
Stokes correction factor 
= empirical parameter for subcritical Reynolds numbers 
= empirical parameter for Mach effects at nexus 
= empirical parameter for rarefaction effects 
n = Knudsen number 
= length scale 
= Mach number 
= particle mass 
= number of particles in a group 
Reynolds number 
= radial coordinate 
= molecular speed ratio of the gas 
gas velocity around the particle surface 
= gas velocity not influenced by the particle 
= particle velocity 
= particle velocity relative to the gas 
= axial coordinate 
= specific heat ratio of gas 
= dynamic viscosity of gas 
= density 
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Subscripts 

accom = surface accommodation for noncontinuum effects 
crit = critical value associated with drag crisis 

D = drag 

fm = free-molecular limit 

g = gas characteristic 

Kn = related to Knudsen number 

m-m = mean free path between molecular collisions 
O = incompressible continuum limit 

p = particle characteristic 

Re = related to Reynolds number 

wall = at solid surface 

0 = tangential momentum 

@p = interpolated to particle centroid 


I. Introduction 


WIDE range of aerospace applications includes spheres and 

spherical particles traveling at high speeds and/or in noncontin- 
uum conditions. For example, hypersonic vehicles exposed to rain, 
ice crystals, and other meteorological particles can undergo signifi- 
cant surface damage due to the high-speed impact of these particles 
[1]. Rocket combustion chambers and flows with dusty detonations 
can also have particles at a high Mach number and/or a high Knudsen 
number [2,3]. Submicron particles can appear in supersonic boun- 
dary layers in various airbreathing propulsion systems [4]. Plasma 
sprays also often include supersonic flow regions with small melting 
metal particles for which two-way coupling effects can be important 
[5]. In such systems, the description of the drag is critical in terms for 
both particle trajectory and momentum coupling. Furthermore, a 
spherical shape is often a reasonable assumption. As such, there is 
significant interest to describe the drag coefficient of a sphere in high- 
speed relative flow conditions [6-12]. 

The quasi-steady drag force Fp of a sphere occurs in the opposite 
direction of the relative particle velocity w, which is the difference 
between the particle velocity v and the gas velocity hypothetically 
interpolated to the particle center u@,. The influence of the relative 
velocity magnitude (i.e., |w|) on the drag force for incompressible 
continuum flow depends on the particle Reynolds number Re,. 
For supersonic and hypersonic flows past a particle, the additional 
effects of compressibility are characterized by the particle Mach 
number M_,, whereas the effects of rarefaction are characterized by 
the particle Knudsen number Kn,. These three nondimensional 
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parameters are related to the particle diameter d, the gas density pg, 
the gas viscosity 4p, the speed of sound a,, the gas temperature T,, 
the gas specific heat ratio y, the gas constant R,, and the gas mean 
free path between molecular collisions J,,_,,, as follows: 


W =V-le, (la) 
d 
Re, = AAE (1b) 
Hg 

m =l ll to 

” ag YR eT. 

i my [M 
Kn, = "s A ld 
Ks d (Re) 


The quasi-steady drag force Fp magnitude can also be expressed 
nondimensionally in terms of a drag coefficient Cp or in terms of the 
Stokes drag correction f as 


1 
Fp = -g 7Pglw|wd* Cp = -3ndu, f w (2a) 
Co= f (2b) 
aan Re, 


The Stokes correction factor f is defined as the ratio of the drag 
coefficient at any Reynolds number relative to that for the creeping 
flow conditions for incompressible continuum flow for which 
Ch = 24/Re,, i.e, f = 1 when Re, <0.01 and Kn, <0.01 
[1,2,10]. In general, the drag coefficient can also be influenced by 
heat and mass transfer rates between the particle and the gas as well as 
the effects of particle shape, lift, added mass, and history forces 
(Refs. [13—20]). However, the present study focuses on nonspinning 
spherical particles in a uniform gas with quasi-steady adiabatic 
conditions. In particular, Cp (or f) are investigated in terms of 
Re,, M,, and Kn, [where specifying two of these parameters 
specifies the third parameter according to Eq. (1d)]. The rarefaction 
and compressibility flow regimes can be characterized by the Mach 
and Knudsen numbers as indicated Tables 1 and 2, where the quanti- 
tative values are based on comprehensive review of available theory, 
data and models for the quasi-steady drag coefficient for spheres in a 
variety of conditions [6—12,21—53]. 

Of the three key parameters to be investigated, a Reynolds number 
has the strongest influence on the drag coefficient. The effects of Re, 
for incompressible and continuum flow conditions (M, — 0 and 
Kn, > 0) on flow regimes and Cp are described in Table 1 and 
Fig. 1 [45,46], respectively. The most viscous condition (Re, — 0) is 
often referred to as the Stokes flow condition, where viscous and 
pressure stresses dominate the flowfield (so inertial effects can be 
neglected). The flowfield in this regime has a theoretical solution 
given by f = 1 per Eq. (2b). This is often referred to as the Stokes 
drag regime and results in a linear dependence of drag on the relative 
velocity in Eq. (2a). As the Reynolds number increases, the flow 
around a particle stays laminar, but convective effects lead to a 
reduced region of viscous effects, resulting in a boundary layer 
around a sphere once Re, > 1. This boundary layer is generally 
attached for Re, < 22, and then it has a steady separated region up 


Table 1  Rarefaction flow regimes 


Knudsen range Flow physics 


Kn, < 0.01 Continuum flow 
0.01 < Kn, < 0.1 Slip flow 

0.1 < Kn, < 10 Transition flow 
Kn, = 10 Free-molecular flow 


Table 2 Compressibility flow regimes for 


a sphere 
Mach range Flow physics 
M, < 9.1 Incompressible flow 


0.1 <M, < 0.65 Compressible subsonic flow 
0.65 <M, < 1.2 
25M <5 

My, 25 Hypersonic flow 


Transonic flow 


Supersonic flow 


to Re, < 130. Figure 1 shows that the drag coefficient in this regime 
is larger than that predicted by Stokes drag (Cp = 24/Re,, i.e., 
f = 1). At still higher Reynolds numbers, the laminar wake becomes 
unsteady and transitional (up to Re, < 1500), and then fully turbu- 
lent (Re, > 1500). 

It can be seen that the drag coefficient is nearly constant for 
1500 < Re Js 250,000. This is often described as the Newton regime 
corresponding to conditions for which the boundary layer on the 
sphere stays fully laminar until it separates, leading to a fully turbu- 
lent wake downstream. However, further increases in the Reynolds 
number are known to lead to a sudden drop in the drag coefficient. 
This is a result of the boundary layer transitioning to turbulent while it 
is still attached. This transition results in a fuller boundary layer that is 
more resistant to separation so that the separation point is moved 
downstream. As result, the turbulent wake zone is smaller such that 
the net pressure drag is reduced leading to the well-known “drag 
crisis.’ The Reynolds number when the attached boundary layer 
changes from laminar to turbulent has been called the “critical” 
Reynolds number Re, crit and it can be defined as the Reynolds 
number when Cp first drops below 0.3 for incompressible continuum 
conditions. As such, the condition whereby the attached boundary 
layer on the sphere stays laminar up to separation is termed “‘sub- 
critical” (Re 250,000 for incompressible flow), whereas the 
condition where the boundary layer becomes transitional or turbu- 
lent before flow separation is termed the “supercritical” regime 
[111,1618]. 

Measurements in the subcritical range have shown that Cp is 
approximately independent of Re, (Fig. 1). This nearly constant 
value is called the “critical drag coefficient” Cp crit with a value of 
approximately 0.42. For the modeling of the drag coefficient from the 
Stokes regime to the critical Reynolds number, one may use the fit 
given by Clift and Gauvin [6] for M, > 0 and Kn, —> 0: 


24 0.42 
Cpo =|——CU+ 0.15Re3)  —— 
ma Le. P 1 + (42, 500/Re!,'*) 
for Re, <2 x 10° and Kn, < 0.01 (3) 


Figure | shows that this widely used incompressible continuum 
model for drag is robust and accurate in this regime (to within 6% of 
experiments). In the following, this Clift-Gavin model is used as the 
basis for adding compressible and rarefaction flow effects for con- 
ditions up until the drag crisis. 

When the Knudsen number increases as in Table 2, the flow around 
the particle can no longer be considered continuum. In this case, the 
molecules do not have a sufficiently high collision rate with the 
particle and each other, leading to a difference between the mean 
molecular velocity and the mean particle surface velocity, i.e., a 
partial-slip condition. This phenomenon is referred “accommoda- 
tion,” and the velocity difference results in a decrease in the drag 
relative to the continuum conditions. For slip flow, the noncontinuum 
effect is weak and can be molded with the Cunningham slip boundary 
condition [21]. Knudsen numbers much greater than unity are gen- 
erally considered “‘free-molecular” flow conditions where the mole- 
cules impact and reflect from the particle surface, not with each other 
(theoretical analysis is again possible in this regime). The “transition” 
regime corresponds to Kn, of order unity and is the most complex 
to theoretically analyze the Boltzmann equation; but, it can be 
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Fig.1 Drag coefficient for a smooth solid sphere at various Reynolds numbers for incompressible continuum flow showing the model of Eq. (3) compared 


with experimental data reported by White [45] and Hoerner [46]. 


addressed numerically with the direct simulation Monte Carlo 
(DSMC) method [22,23]. 

Compressible flow features arise that modify the flowfield 
around the particle once the Mach number increases and the flow 
leads to transonic or supersonic regions, with observed examples of 
key flow features shown in Fig. 2 [47]. If one considers subcritical 
conditions for which a thin boundary layer is formed (e.g., 100 < 
Re, < 1,000,000), the flow will reach a sonic speed once the free- 
stream Mach number is approximately 0.65. The sonic location first 
arises just above the boundary layer at an azimuthal position 0 = 
90 deg away from the leading-edge stagnation point, which is 
consistent with potential flow theory. At still higher freestream 
speeds, the flow becomes transonic with the appearance of a lambda 
shock wave on the particle surface just upstream of 0 = 90 deg 
(Fig. 2a), with the flow becoming turbulent and separated at the 
leading leg of the lambda shock. Once the freestream flow is 
supersonic (M, > 1), the lambda shock is replaced by a standoff 
normal shock in front of the sphere, which transits into an oblique 
shock further downstream (Fig. 2b). The flow separates at just 
upstream of 0 = 90 deg and leads to a weak recompression shock 
at this location. At still higher Mach numbers where the flow is 
hypersonic, the standoff shock becomes very close to the particle; it 
begins to interact with the boundary layer and to create a shock layer 
(Fig. 2c). Further downstream, and oblique shock with a very 
shallow angle forms while the flow separation shear layer becomes 
steadier (because it is more controlled by gas dynamics than turbu- 
lent instabilities). The result of this compressibility is a higher 
contribution of the pressure drag generally leading to a higher drag 
coefficient. In particular, high subcritical Reynolds numbers (e.g., 
10, 000 < Re, < 1,000,000) result in a rapid rise of Cp in the 
transonic regime with a maximum at M ~ 1.4, followed by a 
general reduction in the remaining supersonic regime of 
(1.4 <M, <5) and a relatively constant Cp for hypersonic con- 
ditions (5 < M, < 10) that is still higher than that for incompress- 
ible flow [24]. 

The increase in a Mach number (compressibility) usually leads to 
an increase in the drag coefficient, whereas the increase in a Knudsen 
number (rarefaction) leads to a decrease in the drag coefficient. In 
particular, Loth [24] proposed that compressibility effects dominate 
for Re, > 45 (whereby Cp increases with increasing M ,), whereas 
rarefaction effects dominate for Re, < 45 (whereby Cp decreases 
with increasing M „). In between, it was suggested that the competing 
effects cancel each other out, leading to a “nexus” condition, whereby 
a drag coefficient is independent of a Mach number (and a Knudsen 
number) for Re = 45, i.e., 





Fig. 2 Shadowgraphs and schlieren photographs [47] showing flow 
around spheres at subcritical Reynolds numbers for a) transonic conditions 
(M, ~ 9.8), b) supersonic conditions (M, ~ 2.5), and c) hypersonic condi- 
tions (M, ~ 7.6). 


Cp ~ 1.63 for Re, =45 atallM, and Kn, (4) 


Based on this nexus hypothesis, a rarefied model was developed 
for Re < 45 andacompressible model for Re > 45. Extensive review 
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of the data shows that this model was found to be superior in accuracy 
and robustness as compared to previous models. 

Subsequent to the Loth model [24], Parmar et al. [25] developed a 
new model based on the dataset of Bailey and Hiatt [8] for transonic 
and some supersonic conditions, which provided improved accuracy 
when compared to this particular experimental dataset. However, the 
Parmar et al. model is limited to M, < 1.75 and Kn, < 0.01 (no 
hypersonic conditions) and employs a nonlinear interpolation with 
50 empirically determined coefficients. Furthermore, it was not 
evaluated against much of the previous experimental data reviewed 
for the Loth model. 

Based on the preceding, there is a motivation to develop a drag 
model for spheres that is robust for a large range of Mach numbers 
and should consider all available experimental data, including more 
recent experiments and simulations of particle drag in the supersonic 
and rarefied conditions. Ideally, such a model would be based 
on a modest number of empirical coefficients and would provide 
improved accuracy. Furthermore, such a model should consider 
combinations of Mach and Knudsen numbers where no previous 
data exist. In particular, there are data in the nexus region of Eq. (4) 
based on experimental data with substantial scatter, and there have 
not been high-fidelity simulations of the drag near the nexus for the 
verification of the accuracy of the nexus assumption. Therefore, there 
is aclear need for additional data into this region to better understand 
the trends where effects of both compressibility and rarefaction 
strongly compete. 

The objectives of the current study are threefold. The first objective 
is to conduct a comprehensive data review that not only includes the 
previous experimental data used by Loth [24] and Parmar et al. [25] 
but also incorporates more recent experimental and computational 
data that have since become available. A second objective is to 
conduct new simulations to add additional detail to more accurately 
describe the combined impact of compressibility and rarefaction in 
the vicinity of the hypothetical nexus condition. The third objective is 
to evaluate the previous drag models relative to a full comprehensive 
review of previous data and new simulation data in order to assess 
their accuracy and robustness, as well as to allow the development of 
a new Set of models that could be even more robust and accurate. To 
the authors’ knowledge, this is the first such published study to 
address these three objectives. 


II. Methodology 


A. Comprehensive Data Review 


A detailed survey of all experimental data was conducted. This 
included the data assembled by Loth [24] and Parmar et al. [25], 
resulting in the first study to combine these two experimental datasets. 
However, a key objective of the present study was to compile all 
available high-fidelity data obtained by experiments and simulations. 
Therefore, measurements of high-speed drag coefficients were 
included from several additional sources, including Jacobs [26], Roos 
and Willmarth [27], Short [28], Spearman and Braswell [29], and 
Theofanous et al. [30]. As a result, the experimental data represent 
the most comprehensive review of supersonic and hypersonic sphere 
drag data. In addition, recent direct numerical simulations of con- 
tinuum flow by Nagata et al. [31,32] were also included. However, 
there is a scarcity of high-confidence data in the vicinity of the 
proposed nexus (Re, of around 45 and M, > 0.6). Therefore, addi- 
tional simulations for flow over a sphere were conducted using the 
methods described in the following, and drag coefficients were 
obtained in the regime where both rarefaction and compressibility 
are important. 


B. DSMC Computational Approach 


DSMC was used to provide the drag for conditions with significant 
rarefaction (Kn > 0.1). The simulations were completed using an 
open-source DSMC code developed by Sandia National Laborato- 
ries: SPARTA [33]. Herein, only laminar steady flows are considered, 
which limits the Reynolds number to be less than 45 and allows 
axisymmetric solutions. 


The axisymmetric domain was defined by a radial coordinate r and 
a streamwise coordinate x, as shown in Fig. 3a. An inflow boundary 
was set at 10 diameters upstream (x = —10d), whereas specular 
reflection and outflow boundaries were set at the outer radial boun- 
dary (r = 10d) and at the downstream boundary (x = 10d), respec- 
tively. A periodic boundary condition was imposed in the azimuthal 
direction, consistent with axisymmetric boundary conditions. The 
inflow temperature was set at 293 K, and the isothermal boundary 
condition with 293 K was imposed on the surface of the sphere was 
set as at 293 K. For velocity, a Maxwell boundary condition was 
applied on the surface of the sphere. The accommodation coefficient 
Caccom 1S the fraction of molecules that undergo a diffuse reflection 
upon impacting the particle surface, whereas the remaining fraction 
undergoes specular reflection. In the DSMC predictions, Caccom = 
0.9: a value that is consistent with a smooth surface and experimental 
results [24]. 

The baseline grid used two different resolutions, as shown in 
Fig. 3b, where the near-field grid spacing in all directions was one- 
fourth that of the far-field grid spacing. The grid width and time step 
were Set to be less than the mean free path and the mean collision 
time, respectively. The quantities of the flowfield were sampled as the 
averaged values during 5000 steps to get one averaged flowfield. 
After the flowfield becomes statistically steady, the averaged flow- 
field was computed again using more than 20 averaged flowfields. 
Examples of the velocity and number-density distributions are shown 
in Fig. 4. The drag coefficient was calculated using flow data after the 
flowfield achieved a statistically steady state, and grid and time-step 
resolution studies indicated that the drag values converged to within 
1% for all simulations. 


C. Navier-Stokes Computational Approach 


To consider conditions with weak or no rarefaction, additional 
simulations were conducted with the full Navier-Stokes equations, 
with all viscous effects fully resolved. This approach is often referred 
to as Direct Numerical Simulation (DNS) and can be used to treat 
flows in the laminar, transitional, or turbulent regime. In the 


Specular reflection 


-X 
Inflow boundary 
(Particle emit surface) 


+x 
Outflow boundary 


Surface of sphere 
Maxwell reflection 





d 
r=0 
a) Axi-symmetric boundary 





b) 


Fig. 3 DSMC flow plane for the streamwise direction x and radial 
direction r showing a) boundary conditions for this plane (periodic 
boundary conditions used in the azimuthal direction) and b) computa- 
tional grid showing course and fine grid regions. 
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Fig.4 Example DSMC results for Re, = 45 and Kn, = 0.242 (M, = 1.2) in the streamwise-radial plane showing contour distributions of a) velocity 
relative to the sphere (in meters per second) and b) number density (per cubic meter). 


continuum regime (Kn, < 0.01), these equations can be employed 
with a no-slip condition. For weak rarefaction conditions 
(0.01 < Kn, < 0.1), the surface velocity conditions can be treated 
using a Maxwell slip relationship [34]. In particular, the difference 
between the local wall tangential velocity U way and that of the fluid U 
for a spherical coordinate system is related to the distance above the 
surface (r—r,), the radial velocity gradient, and the tangential 
momentum coefficient cg as 





Cod oU 
Usa = U + = $ — Kn, (5a) 
p 
2— 
Co = Caccom (5b) 
Caccom 


For weak rarefaction, the focus was on laminar steady flows with 
Re, < 200 for subsonic conditions and Re, < 300 for transonic 
conditions to ensure steady axisymmetric conditions. 

Adiabatic conditions were also assumed, which is consistent witha 
particle temperature equal to the surrounding gas temperature. The 
particle temperature can be different from the ambient temperature if 
the particle is undergoing rapid heating or cooling (e.g., due to 
injection conditions, shock passage, etc.). However, the effect of 
temperature difference on drag is generally small for the Reynolds 
numbers simulated based on experimental data [7,8]. 

Discretized governing equations were solved by an in-house com- 
putational code with high-order schemes, boundary conditions, and 
computational grids where the computational code has been vali- 
dated in detail for drag over spheres [31,32]. The convection terms 
were evaluated by the sixth-order adaptive central and upwind 
weighted essentially nonoscillatory scheme WENOCU6-FP of Non- 
omura et al. [35], but the central difference component was replaced 
by the splitting type proposed by Pirozzoli [36] in order to stabilize 
the calculation. The viscous terms were evaluated by the sixth-order 
central difference method. The time integration is conducted by the 
third-order total variation diminishing Runge—Kutta method pro- 
posed by Gottlieb and Shu [37]. 


III. Results 
A. Quasi Nexus 


The previously proposed nexus condition suggested that the drag 
coefficient was independent of a Mach number for Re, = 45 [per 
Eq. (4)]. This trend is consistent with results that show that higher 


Reynolds numbers generally yield an increase in Cp as M , increases 
(due to compressibility effects that create a higher-pressure drag over 
the surface), whereas lower Reynolds numbers generally lead to a 
decrease in Cp as M, increases (due to rarefaction effects that result 
in a slip velocity on the surface). 

To investigate the validity of the nexus, all available data including 
new simulations conducted for this study were collected for a range of 
Mach numbers at Re, ~ 45, as shown in Fig. 5 [48,49]. Although 
there is considerable spread and uncertainty in the experimental data 
(gray and black filled symbols), the computational results (open 
symbols) clearly indicate a transonic bump that is inconsistent with 
the Loth model [24], which is shown as a dashed line. In particular, 
the drag coefficient increases in the subsonic range from the incom- 
pressible value of 1.63 for M, <0.1 to Cp ~2 at M, ~1. For 
supersonic Mach numbers, the drag coefficient decreases and then 
returns approximately to the incompressible values of 1.63 for hyper- 
sonic conditions (M, > 5). These results indicate that the nexus 
proposed in Eq. (4) is not strictly realized. Instead, a more accurate 
description for Re, near 45 is a “quasi nexus” where the influences of 
a Mach number are minimized but still significant. In particular, the 
expected continuum flow increases in Cp due to compression tended 
to be mostly (but not completely) canceled by the decreases in Cp due 
to rarefaction in the vicinity of Re, of 45. In particular, the present 
simulations found that the wave drag effects tend to increase the 
pressure on the fore side of the sphere (thereby increasing drag), 
whereas rarefaction effects reduce the velocity gradient on the surface 
due to partial slip (thereby decreasing drag). This is an important 
aspect in constructing empirical models for the drag coefficient. 

This figure also shows the empirical models for the compressible 
regime and for the rarefied regimes, which are seen to be well bridged 
at this Reynolds number (1.e., within 1% of each other) and reason- 
ably reflect the transonic bump trends described earlier in this paper. 
These two empirical models are discussed in the following two 
sections. However, there is a wide variation in data at this particular 
Reynolds number. Experimental measurements are notoriously hard 
to obtain with high accuracy because the entire flow can be highly 
unstable due to the shock movement interactions with the boundary 
layer. In particular, small surface variations or upstream unsteadiness 
can cause significant shock movement, thereby changing the surface 
pressure over the sphere. In addition, transonic flows are highly 
sensitive to flow blockage effects, with the present simulations 
indicating a boundary condition of 10 diameters (Fig. 3) as being 
critical for this condition, which is more than any other Mach number 
conditions. This may be one of the reasons that there have been no 
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e Bailey and Hiatt [8] (Re = 41-49) 


= Roos and Willmarth [27] (Re = 45-54) 
a Wegener and Ashkenas [45] (Re = 35-46) 


e Aroesty [7] (Re = 35-48) 
a Carlson and Hoglund [49] (Re = 40) 
4 White [45] (Re = 48) 


o Present DNS with slip (Re = 45) 
oO Present DSMC (Re = 45) 
à Nagata and Ashkenas [48] DNS (Re = 50) 
Present compressible model 
Present rarefied model 
= = = Loth model [24] 





simulations [31], present DNS and DSMC simulations, nexus model [24] [Eq. (4)], as well as present models for the compressible flow regime 


[Eqs. (6-8)] and the rarefied regime [Eqs. (9-15)]. 


previous detailed simulations for this condition before the present 
work. As such, further work is needed to improve the understanding 
and modeling of the quasi nexus. 


B. Compression-Dominated Regime 


The influence of a Mach number on the drag coefficient is well 
known at high Reynolds numbers in the subcritical regime, as dis- 
cussed in the Introduction (Sec. I). As expected, the compression- 
dominated regime generally yields increased drag as the Mach 
number increases, which then settles to a near-constant value. Based 
on Refs. [24] and [25], the effects of compressibility at high Reynolds 
numbers (e.g., 10,000 < Re ps 1,000,000) can be expressed in 
terms of a drag coefficient ratio defined as 





Cp crit Cp crit 
Cy = =n (6) 
" Cperixm—o 0.42 


This expression uses a baseline (incompressible) critical drag 
coefficient of 0.42, which is consistent with the value at high 
subcritical Reynolds numbers as shown in Fig. 1. However, the 
current study sought to review all available experimental data and 





to accurately describe the influence of M, on Cy. This updated 
collection of data is shown in Fig. 6. It can be seen that the models 
of Loth [24] and Parmar et al. [25] reasonably captured the trends in 
Cy, but the Parmar et al. model tends to overpredict at lower Mach 
numbers. It should be noted that the previous models were only based 
on the data reported by Hoerner [46]. Based on this updated collec- 
tion of all data available (which differs slightly), a new model is 
proposed for Cy as 


Cy = 1.65 +0.65tanh(4M, —3.4) for M,<1.5 (7a) 


Cy = 2.18 — 0.13 tanh(0.9M, —2.7) forM,>1.5 (Tb) 


This new model is simpler than that proposed by Loth [24] and is the 
most accurate of the models shown in Fig. 6. 

To investigate the influence of Mach number between the quasi- 
nexus point (Fig. 5) and the high-Re, regime (Fig. 6), data were 
collected (including present simulations) as a function of a Reynolds 
number for a set of Mach numbers as shown in Fig. 7 [46-52]. For the 


nearly incompressible case of M, ~ 0.2 (Fig. 7a), there is the local 
minimum of Cp in the range of 1000 to 10,000 beyond which the drag 


Roos and Willmarth [27] 
= Bailey and Hiatt [8] 

Short [28] 

White [45] 

Bailey and Starr [38] 


Spearman and Braswell [29] 
Hoemer |46] 

Present model 

". Parmar et al. [25] 
Loth model [24] 
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Bailey and Hiatt [51] (M = 0.23) 
Goin and Lawrence [52](M= 0.20 & 0.33) 
Theofanous et al. [30] (M ~ 0.28) 
Bailey and Starr [38] (M ~ 0.28) 
Roos and Willmarth [27] (M = 0.25) 
Jacobs [26] (M = 0.2) 
Present DNS with slip (M = 0.20) 
Present model 
-= = - Loth model [24] 
Farmar et al. [25] 


Bailey and Hiatt [51] (M = 0.49) 
Goin and Lawrence [52] (M = 0.60) 
Theofanous et al. [30] (M ~ 0.52) 
Fresent DNS (M = 0.50) 
Present DNS with slip (M = 0.50) 
Fresent model 
Farmar et al. [25] 

=- = -Loth model [24] 


Bailey and Hiatt [51] (M = 0.75) 
Goin and Lawrence [52] (M = 0.75) 
Short [28] (M = 0.75) 
Theofanous et al. [30] (M ~ 0.72) 
Present DNS with slip (M = 0.80) 
Nagata et al. [31] DNS (M = 0.80) 
Present model 

rrenen - Parmar et al. [25] 

= = -Loth model [24] 





Fig. 7 Drag coefficient in compressible regime (Re, > 45) as a function of Reynolds number based on experimental and computational data 
[7,8,26—32,38,46—52] as well as present simulations and empirical models including the present compressible flow model [Eqs. (6-8)] for a) M, ~ 0.25, 
b) M, ~ 0.5, c) M, ~ 0.75, d) M, ~ 0.83, e) M, ~1,f) Mp ~ 1.5, g) M, ~ 2, h) M, ~4,i) M, ~5.5, and j) M, ~ 10. 
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Fig. 7 (Continued) 
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Fig.7 (Continued) 


rises to its critical value for the Newton regime of 0.42 (as seen in To describe the compressible flow trends in the drag coefficient 
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Fig. 1). Ateven higher Reynolds numbers, the drag drops again due to 
transition of the attached boundary layer from laminar to turbulent. 
The associated drag crisis is shown to occur at Re, crit ~ 150,000 in 
this case. As the Mach number increases to transonic conditions 
(Figs. 7b—7e), the drag coefficient is elevated at all Reynolds num- 
bers. In addition, the local minimum of the Cp in the range of 1000 to 
10,000 is eliminated and the drag crisis is delayed at higher Reynolds 
numbers (Re, crit > 800,000 for M, > 0.75). As the Mach number 
becomes supersonic and hypersonic (Figs. 7f{—7]), the drag coeffi- 
cient is seen to continually reduce as the Reynolds number increases. 
In addition, the drag crisis continues to be delayed (e.g., Rep crit > 
1,500,000 for M, > 2) based on experimental results [38]. This is 
attributed to the increasing stability of the boundary layer on the 
sphere surface as the Mach number rises, as has been demonstrated 
experimentally and theoretically for flat-plate boundary layers 
[39,40]. Although there are insufficient experimental data to accu- 
rately provide a model of Re, crit for various M,, this influence 
should be kept in mind when employing high-speed drag coefficient 
models. 





based on all the available experiments and simulations with 
Re, > 45, a modified Clift-Gauvin drag expression was developed 
with Eq. (7) as 


24 
Cp = Re. (1 + 0.15Re°8’)Hy 
P 
n 0.42Cy 
1 + (42, 500/Rep ®™) + (Gy /Re®) 
for 45 < Re, < Rep crit (8a) 


Gy = 166M3 + 3.29M} — 10.9M, +20 for M,<0.8 (8b) 
Gy =5+40M;3 for M, > 0.8 (8c) 


Hy = 0.0239M}, + 0.212M?, —0.074M, +1 for M, <1 (8d) 


è Millikan [12] (Oil Droplets) 
à Millikan [12] (Mercury Droplets) 
m Millikan [12] (Shellac Droplets) 
© Benson etal. [42] DSMC 

Clift and Gauvin [6] fit 


===- Phillips theory [9] 


Fig.8 Stokes drag correction as a function of Knudsen number for experimental and computational data [12,42] along with the Clift and Gauvin model of 


Eq. (10). 
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l 
Se E for M, > 1 (8e) 


Figures 7a—7j illustrate that the present model is more accurate 
and robust for supersonic and hypersonic conditions than the 
previous models of Loth [24] and Parmar et al. [25]. Furthermore, 
the present formulations of | Gy and Hy are simpler i in form than 
those previously proposed by Loth [24]. The quantitative dif- 
ferences between the various model predictions compared to the 
data from experiments and high-fidelity numerical simulations are 
discussed later. 


C. Rarefaction-Dominated Regime 


The effect of rarefaction becomes important for finite Knudsen 
numbers, and this effect tends to dominate the drag coefficient trends 
of Re, < 45. When the flow is in continuum (Kn, — 0), the drag 
force for a spherical particle given by Eq. (3) is reasonable for 
Re, <45. For finite Knudsen numbers in creeping flow 
(Re, > 0), the Cunningham correction factor is defined as 


Fp Re-0 
= 9 


For a small Mach number (M, « 1), Clift and Gauvin [6] developed 
a simple empirical model that is consistent with cg¥1.22 and the 
approximate theory of Phillips [9] (which bridged the theoretical 
result by Basset [10] for Kn, « 1 and that of Epstein [11] for Kn, 1) as 


l 


ee 10 
fxn 1 + Kn [2.514 + 0.8 exp(—0.55/Kn,)] a 


As shown in Fig. 8, this model agrees with the theoretical 
results and experiments of Milliken [12,41] and the DSMC of 
Benson et al. [42] that showed a decrease in drag with increasing 
Kn,. To extend Eq. (10) to include finite Reynolds numbers in the 
rarefaction-dominated regime (Re, < 45), the Knudsen number 
effect can be modified to include a Schiller-Naumann correction 
[24] as 


4 
Cp.Kn.Re = a) + 0.15Re}°") fen (11) 
P 


It should be noted that Eq. (11) is limited to weak compress- 
ibility effects (M,.«1). 

To incorporate rarefied conditions with high-Mach-number con- 
dition effects, one may instead consider the free-molecular limit, 
which can be characterized by the molecular speed ratio s defined 
by the Mach number and specific heat ratio as 


s=M,Vy/2 (12) 


Assuming equal tangential and normal accommodation coeffi- 
cients and equal particle and gas temperatures, the free-molecular 
compressibility limit for s > 1 as obtained by Stalder and Zuric [43] 
and Schaff and Chambré [44] is 


(1 + 2s”) exp(—s? im (4st + 4s? — 1)erf(s) 
Cp tim = > + HMMM 13 
ae S/T 254 gr = Ja os 
In the limit of s > 00, Cp fm > 2, which can be taken as the lower 
limit of the rarefied drag coefficient for Re p<i Herein, this free- 
molecular limit can be empirically corrected for finite Reynolds 
numbers as 


Cp, fm 


~ 14 [(Cppm/Iu) — 1],/Re,/45) 


CD, fm,Re (14a) 


(14b) 





2,20 = it for M, <1 
ag 1.6 + 925 + O1 4 OM for M, > 1 


Compared to those previously proposed by Loth [24], the present 
forms of Cp fm.re and J y are simplified via the assumption of equal 
gas and particle temperatures but also fix a typo for Cp fm re that 
should have included the last specular reflection term in Eq. (13). 


© Present DNS with Slip (M = 0.8) 
- M=0, Kn=0 (Clift and Gauvin [6]) 
Present model 
Loth model [24] 


Aroesty (M -7) 

Carlson and Hoglund [49] (M = 2) 
Bailey and Hiatt [8] (M ~ 2) 
Present DSMC (M = 2) 

Present DNS with Slip (M = 2) 
Present model 

Loth model [24] 


() 
0.00] 


4 


Aerosty [7] (M ~ 3) 
Macrossan [53] DSMC (M = 3) 
Present DNS with Slip (M = 3) 
Present DSMC (M = 3) 

Present model 

Loth model [24] 


® Aerosty [7] (M~4) 
O Present DNS with Slip (M = 4) 





Fig. 9 Drag coefficient in rarefaction regime (Re, < 45) as a function 
of Reynolds number based on previous data [7, ‘8, 49,53], as well as 
present simulations and empirical models [Eqs. (9-15)] for a) M, ~ 0.8, 
b) M, ~2, ¢) M, ~ 3, d) M, ~4, e) M, ~5, and f) M, ~ 10. 
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à Macrossan [53] DSMC (M = 5) 
© Present DSMC (M =5) 
Present model 


| A Macrossan [53] DSMC M = 10) 
© Present DSMC (M = 10) 
Present model 





Fig.9 (Continued) 


Finally, the weak Mach-number-limit of Eq. (11) can be blended with 
the high-Mach-number limit of Eq. (14) in the rarefaction-dominated 
regime using the empirical expression of Loth [24]: 


= Cp. Kn.Re M$ CD fm.Re 
1+M$ 1+M$ 





D for Re, < 45 (15) 


The preceding set of equations integrates the combined influence 
of the Reynolds, Mach, and Knudsen numbers for the rarefaction- 
dominated regime and are simpler in form than the model of 
Loth [24]. 

Figures 9 and 10 show the overall behavior for the rarefaction 
regime with an increasing Mach number. In general, increasing 
rarefaction results in a decreasing drag coefficient, which is consis- 
tent with increasing the slip conditions on the particle surface. In 
addition, it can be seen that the drag coefficient tends to a constant 
value for Re, « 1 for all the supersonic and hypersonic conditions that 
were evaluated (Figs. 9a—9g) [53]. For the hypersonic conditions, this 
tendency yields Cp > 2 as Re, — 0, which is consistent with the 
free-molecular limit (shown by the dashed lines in Fig. 10). 

In general, the present empirical model given by Eqs. (10-15) 
predicts all the trends observed by the experiments and simulations, 
including the new simulations conducted for this study. The model 
predictions are effectively the same for hypersonic flow (M, > 5) as 
those of Loth [24] but are more accurate for subsonic and supersonic 
flows. In general, the model predicts the limiting value of a drag 
coefficient for Re, « 1 observed in the data for supersonic and hyper- 
sonic conditions (Fig. 10). However, it should be noted that the model 
also predicts a limiting behavior for M, = 0.8, as shown in Fig. 9a 
(Cp > 6.2 as Re, — 0), but there are no experiments or simulations 
that confirm this result. As a result, it is recommended that simula- 
tions and experiments should be conducted to better characterize the 
influence of Mach number on the drag coefficient in the subsonic 
limit of creeping flow at rarefied conditions. However, the current 
models accurately and robustly predict all the trends and values given 
by the available data in this comprehensive review. The overall trends 
are Shown in Fig. 11, including the quais-nexus. 


D. Accuracy Assessment of Models and Overall Drag Behavior 
Figures 6-10 show that the present model for the drag coefficient of 

supersonic and hypersonic flows past spheres is generally accurate 

and robust in terms of predictive performance for a wide range of 


o M=2 DSMC (Macrossan [53]) 
O M=3 DSMC (Macrossan [53]) 
å M=5 DSMC(Macrossan [53]) 
© M=10 DSMC (Macrossan [53]) 
© M=20 DSMC (Macrossan [53]) 


—— M=0, Kn=0 (Clift and Gauvin [6]) 
Rep<<1 and Knp<<1 theory: CD',fm 
"| ——= M>>1(s>>1) fit: CD,fm,Re 
Present fit: CD fm.M_Re 





Fig. 10 Drag coefficient (Cp) variations in rarefaction regime 
(Re, < 45) showing limits for incompressible continuum flow (thick solid 
line), free-molecular flow (long dash line), and hypersonic creeping flow 
(short dashed line), along with typical DSMC results [53]. 


Mach, Knudsen, and Reynolds numbers. To quantitatively assess the 
accuracy of the present model and those of Loth [24] and Parmar et al. 
[25], a detailed error analysis relative to all the experimental and 
computational results presented in this paper was completed. Since 
the error was based on the difference between each data point pre- 
sented, the error included uncertainty associated with the data itself. 
For example, the significant variation in the experimental data at 
transonic conditions (e.g., Fig. 6) meant that there was a lower bound 
on the accuracy of any empirical model that could be proposed. 

The net errors for three Mach regimes are shown in Table 3, where 
it can be seen that the present model is approximately twice as 
accurate as either the Loth model [24] or the Parmar et al. model 
[25] for moderate Mach numbers (M, <2) in compressible- 
dominated flow and is generally more accurate than the Loth model 
[24] for both higher Mach numbers (M, > 2) and for rarefied con- 
ditions (Re, < 45). Since those two previous models were found 
to be superior to others presented previously [39,40], the present 
model appears to be the most accurate and robust empirical model 
yet developed for compressible-dominated and rarefied-dominated 
conditions. 


IV. Conclusions 


A comprehensive review of all relevant experimental data 
was completed to characterize the influence of relative Mach, 
Knudsen, and Reynolds numbers on the drag coefficient of a sphere 


Table3 Average root-mean-square difference for drag models based 
on all available data? 


Moderate Mach High Mach 
Rarefied flow (Re, >45,M <2), (Re, 2 45, M > 2), 
Model (Re, < 45), % % % 
Loth [24] 10.7 9.9 3.6 
Parmar [25] N/A 10.0 N/A 
Present model 6.1 4.9 3.4 


“N/A denotes “not applicable.” 
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Fig.11 Rarefaction and compression effects on drag of a spherical particle based on updated models for compressibility and rarefaction, showing effect 


of transonic bump on quasi nexus. 


at supersonic and hypersonic conditions. At high Reynolds numbers, 
the results indicated that the drag crisis was delayed to higher Reynolds 
numbers, indicating Re, crit > 10° for supersonic conditions. 

In addition, new Navier-Stokes and DSMC simulations were 
conducted to investigate the behavior of the drag coefficient at 
Reynolds numbers near 45. The overall combination of the super- 
sonic and hypersonic regimes provided by the data (and predicted by 
the present model) are shown in Fig. 11, where the compression- 
dominated region (Re, < 30) yields an increase in Cp as M, 
increases; whereas the rarefaction-dominated (Re, > 60) yields a 
decrease in Cp as M, increases. In between these regimes, there is a 
quasi nexus with a weak transonic bump occurring that bridges the 
rarefaction-dominated regime and the compression-dominated 
regime. The present simulations are the first to quantitatively describe 
this transonic bump for this Reynolds number, but the wide spread of 
experimental data indicates further work is needed to better under- 
stand and model this quasi nexus. 

Empirical models were also developed, taking advantage of the 
comprehensive experimental dataset, new simulation results, and an 
improved understanding of the quasi-nexus behavior, especially for 
transonic conditions in the compression-dominated region. The new 
empirical models are as much as twice more accurate than previous 
models and served to provide an overall description of the drag 
coefficient for a wide range of Mach, Knudsen, and Reynolds num- 
bers. However, additional simulations and experiments are recom- 
mended to better understand the fluid physics and the drag coefficient 
in the quasi-nexus regime, as well as in the influence of Mach number 
in the subsonic limit creeping flow for rarefied conditions. 
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